home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / specfunc / bessel_Knu.c < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  4.7 KB  |  165 lines
  1. /* specfunc/bessel_Knu.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* Author:  G. Jungman */
  21.  
  22. #include <config.h>
  23. #include <gsl/gsl_math.h>
  24. #include <gsl/gsl_errno.h>
  25. #include <gsl/gsl_sf_exp.h>
  26. #include <gsl/gsl_sf_gamma.h>
  27. #include <gsl/gsl_sf_bessel.h>
  28.  
  29. #include "error.h"
  30.  
  31. #include "bessel.h"
  32. #include "bessel_temme.h"
  33.  
  34. /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
  35.  
  36. int
  37. gsl_sf_bessel_Knu_scaled_e(const double nu, const double x, gsl_sf_result * result)
  38. {
  39.   /* CHECK_POINTER(result) */
  40.  
  41.   if(x <= 0.0 || nu < 0.0) {
  42.     DOMAIN_ERROR(result);
  43.   }
  44.   else {
  45.     int N = (int)(nu + 0.5);
  46.     double mu = nu - N;      /* -1/2 <= mu <= 1/2 */
  47.     double K_mu, K_mup1, Kp_mu;
  48.     double K_nu, K_nup1, K_num1;
  49.     int n;
  50.  
  51.     if(x < 2.0) {
  52.       gsl_sf_bessel_K_scaled_temme(mu, x, &K_mu, &K_mup1, &Kp_mu);
  53.     }
  54.     else {
  55.       gsl_sf_bessel_K_scaled_steed_temme_CF2(mu, x, &K_mu, &K_mup1, &Kp_mu);
  56.     }
  57.  
  58.     /* recurse forward to obtain K_num1, K_nu */
  59.     K_nu   = K_mu;
  60.     K_nup1 = K_mup1;
  61.  
  62.     for(n=0; n<N; n++) {
  63.       K_num1 = K_nu;
  64.       K_nu   = K_nup1;
  65.       K_nup1 = 2.0*(mu+n+1)/x * K_nu + K_num1;
  66.     }
  67.  
  68.     result->val = K_nu;
  69.     result->err = 2.0 * GSL_DBL_EPSILON * (N + 4.0) * fabs(result->val);
  70.     return GSL_SUCCESS;
  71.   }
  72. }
  73.  
  74.  
  75. int
  76. gsl_sf_bessel_Knu_e(const double nu, const double x, gsl_sf_result * result)
  77. {
  78.   gsl_sf_result b;
  79.   int stat_K = gsl_sf_bessel_Knu_scaled_e(nu, x, &b);
  80.   int stat_e = gsl_sf_exp_mult_err_e(-x, 0.0, b.val, b.err, result);
  81.   return GSL_ERROR_SELECT_2(stat_e, stat_K);
  82. }
  83.  
  84.  
  85. int
  86. gsl_sf_bessel_lnKnu_e(const double nu, const double x, gsl_sf_result * result)
  87. {
  88.   /* CHECK_POINTER(result) */
  89.  
  90.   if(x <= 0.0 || nu < 0.0) {
  91.     DOMAIN_ERROR(result);
  92.   }
  93.   else if(nu == 0.0) {
  94.     gsl_sf_result K_scaled;
  95.     /* This cannot underflow, and
  96.      * it will not throw GSL_EDOM
  97.      * since that is already checked.
  98.      */
  99.     gsl_sf_bessel_K0_scaled_e(x, &K_scaled);
  100.     result->val  = -x + log(fabs(K_scaled.val));
  101.     result->err  = GSL_DBL_EPSILON * fabs(x) + fabs(K_scaled.err/K_scaled.val);
  102.     result->err += GSL_DBL_EPSILON * fabs(result->val);
  103.     return GSL_SUCCESS;
  104.   }
  105.   else if(x < 2.0 && nu > 1.0) {
  106.     /* Make use of the inequality
  107.      * Knu(x) <= 1/2 (2/x)^nu Gamma(nu),
  108.      * which follows from the integral representation
  109.      * [Abramowitz+Stegun, 9.6.23 (2)]. With this
  110.      * we decide whether or not there is an overflow
  111.      * problem because x is small.
  112.      */
  113.     double ln_bound;
  114.     gsl_sf_result lg_nu;
  115.     gsl_sf_lngamma_e(nu, &lg_nu);
  116.     ln_bound = -M_LN2 - nu*log(0.5*x) + lg_nu.val;
  117.     if(ln_bound > GSL_LOG_DBL_MAX - 20.0) {
  118.       /* x must be very small or nu very large (or both).
  119.        */
  120.       double xi  = 0.25*x*x;
  121.       double sum = 1.0 - xi/(nu-1.0);
  122.       if(nu > 2.0) sum +=  (xi/(nu-1.0)) * (xi/(nu-2.0));
  123.       result->val  = ln_bound + log(sum);
  124.       result->err  = lg_nu.err;
  125.       result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  126.       return GSL_SUCCESS;
  127.     }
  128.     /* can drop-through here */
  129.   }
  130.  
  131.  
  132.   {
  133.     /* We passed the above tests, so no problem.
  134.      * Evaluate as usual. Note the possible drop-through
  135.      * in the above code!
  136.      */
  137.     gsl_sf_result K_scaled;
  138.     gsl_sf_bessel_Knu_scaled_e(nu, x, &K_scaled);
  139.     result->val  = -x + log(fabs(K_scaled.val));
  140.     result->err  = GSL_DBL_EPSILON * fabs(x) + fabs(K_scaled.err/K_scaled.val);
  141.     result->err += GSL_DBL_EPSILON * fabs(result->val);
  142.     return GSL_SUCCESS;
  143.   }
  144. }
  145.  
  146.  
  147. /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
  148.  
  149. #include "eval.h"
  150.  
  151. double gsl_sf_bessel_Knu_scaled(const double nu, const double x)
  152. {
  153.   EVAL_RESULT(gsl_sf_bessel_Knu_scaled_e(nu, x, &result));
  154. }
  155.  
  156. double gsl_sf_bessel_Knu(const double nu, const double x)
  157. {
  158.   EVAL_RESULT(gsl_sf_bessel_Knu_e(nu, x, &result));
  159. }
  160.  
  161. double gsl_sf_bessel_lnKnu(const double nu, const double x)
  162. {
  163.   EVAL_RESULT(gsl_sf_bessel_lnKnu_e(nu, x, &result));
  164. }
  165.